\chapter{Test Case 1}

\section{Stability and Accuracy of various time integration schemes.}

To assess the stability of various implicit time integration procedures, we have considered a simplified problem in one dimension given below.


\begin{subequations}
\label{Reference_Test_case}
\begin{equation}
\frac{\partial \zeta}{\partial t} +\textbf{U} \dfrac{\partial \zeta}{\partial x} + h \dfrac{\partial^2 \varphi}{\partial x^2} = 0,
\end{equation}
\begin{equation}
\frac{\partial \varphi}{\partial t} + \textbf{U}\dfrac{\partial \varphi}{\partial x} + g\zeta = 0,
\end{equation}
\end{subequations}

Here $\textbf{U}$ represents the mean current velocity and $h$ represents the depth of the ocean. Both $\textbf{U}$ and $h$ have been assumed to be uniform in $x$ direction.

Periodic boundary conditions have been assumed for both $\zeta$ and $\varphi$, and the initial conditions are given as:
\begin{itemize}
 \item $ \zeta  (0,x)$ = cos($\dfrac{4\pi x}{L}$)  
   \item $ \varphi  (0,x)$ = 1
\end{itemize}



Central difference scheme has been used for the spatial discretization. 

\begin{subequations}
\label{test_separate_p}
\begin{equation}
\frac{\partial \vz}{\partial t} + \Szz \vz + \Szp \vp =0
\end{equation}
\label{test_separate_q}
\begin{equation}
\frac{\partial \vp}{\partial t} + \Spz \vz + \Spp \vp = 0
\end{equation}
\end{subequations}

Where the discretization matrices are given by three point stencils as given below.


\begin{align*}
 \Szz  &: \begin{bmatrix} \dfrac{-U}{2 \Delta x} && 0  && \dfrac{U}{2 \Delta x}\end{bmatrix}\\  
 \Szp  &: \begin{bmatrix}  \dfrac{h} {\Delta x^2} && \dfrac{-2 h} {\Delta x^2} && \dfrac{h} {\Delta x^2} \end{bmatrix} \\
 \Spz  &: \begin{bmatrix} 0 && g && 0 \end{bmatrix} \\
 \Spp  &: \begin{bmatrix} \dfrac{-U}{2 \Delta x} && 0  && \dfrac{U}{2 \Delta x}\end{bmatrix}
\end{align*}

The matrices $\Szz$ and $\Spp$ are skew- symmetric, matrix $\Szp$ is symmetric, whereas $\Spz$ is only a diagonal matrix. For the sample problem, following values have been used:

\begin{itemize}
  \item $x-length$ = L =100 m
  \item $\Delta x$ = 0.5 m
  \item U = 1.0 m/sec
  \item h = 50 m
  \item $time-max$ = 100 sec
\end{itemize}

For the combined system of equations represented by \Cref{test_separate_p}, the maximum Eigenvalue obtained is :$\lambda_{max} = 2\sqrt{gh} $. Time step is limited by the CFL condition for the explicit time integration schemes. It is given as:

\begin{equation}
\label{CFL_Number}
  \dt \leq \dfrac{\Delta x}{\lambda_{max}} \leq \dfrac{\Delta x}{2\sqrt{gh}} = 0.0113 \text{ sec} = \dt_{cfl}
\end{equation}

Various time integration schemes as discussed in Chapter 4 have been implemented. In case of Fully Implicit Schemes (Fully Implicit Backward Euler, Fully Implicit Trapezoidal etc.), \Cref{test_separate_p} is solved simultaneously. In case of Semi-Implicit schemes, an iterative procedure is followed to solve \Cref{test_separate_p}.

\subsection{Benchmark Solution using Matlab integrators}
\label{benchmark}

In order to get the benchmark solution, Matlab ODE integrators are used. The system of ODE's given by \Cref{test_separate_p} are integrated by ODE45 and ODE15s functions. ODE45 is an explicit Runge-Kutta method with adaptive time stepping. The time step is modified based on the error estimate provided. ODE15s is a higher order implicit method. For the benchmark, relative error estimate of $1e^{-6}$ is used.

Using ODE15s, we obtain the following solution for the the variable $\zeta$.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15s_benchmark_2}~\\[0.25cm]
  \caption{$\zeta$ variation with length at different times using ODE15s}
  \label{ode15s_benchmark_2}
\end{figure}

The result depends on the Initial condition chosen. If the Initial conditions are smooth, as in \Cref{ode15s_benchmark_2}, the final solution is also smooth.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15s_benchmark}~\\[0.25cm]
  \caption{$\zeta$ variation with length at different times using ODE15s with discontinuity}
  \label{ode15s_benchmark}
\end{figure}

In case if we introduce some discontinuity in the initial condition as shown in in \Cref{ode15s_benchmark} at x=0, wiggles/oscillations are observed with the Central discretization. In order to remove these oscillations, high resolution spatial discretization schemes (flux-limiters or total variation diminishing) are required, which are not considered in the current study.

The comparison of the ODE45 (explicit) and ODE15s (implicit) is shown in \Cref{ode45vs15s}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode45vs15s}~\\[0.25cm]
  \caption{$\zeta$ variation with length at different times using ODE15s and ODE45}
  \label{ode45vs15s}
\end{figure}

In order to compare the two solutions, we define the norm of the difference between the two vectors. Let $x_s$ be the number of grid points in $x$ direction.
\begin{equation}
\label{norm}
 Norm(u,v) = \sqrt{\frac{1}{x_s}\sum_{i=1}^{x_s} (u(i)-v(i))^2}
\end{equation}

For the solutions given in \Cref{ode45vs15s}, the norm as given by \Cref{norm} is $2.0500e^{-4}$. We will now use ODE45 as the benchmark for the explicit time integration solvers given in Chapter 4, and ODE15s as the benchmark for the implicit and semi-implicit time integration methods.


\subsection{Explicit Scheme}

The Leap Frog scheme has been used as the explicit time integration scheme as described in the Chapter 4. For the first time step as the value at two previous time steps in unknown, the Forward Euler method is used. The Leap Frog scheme requires the CFL condition for the stability of the time integration. For time step upto $\dt_{cfl}$ sec, the method is stable. It becomes unstable when the time step is equal or greater than $\dt_{cfl}$ sec.

For various time steps, simulation was run in Matlab for 100 seconds, and the results plotted in \Cref{ode45vsleapfrog}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode45vsleapfrog}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE45 vs Explicit Leap Frog}
  \label{ode45vsleapfrog}
\end{figure}

We can see that as the time step is increased, the accuracy is impacted. In order to compare the results, the norm defined in \Cref{norm} is used and tabulated below.


\begin{table}[H]
  \centering
  \caption{Leap Frog vs ODE45, Error Norm}
  \vspace{0.2cm}
  \begin{tabular}{|c|c|}
  % use packages: color,colortbl
  \hline
  \rowcolor{tcA}
  Time step (sec) & Error Norm\\
  \hline
  0.0113 ($\dt_{cfl}$) & 0.0370\\ \hline
  0.01 & 0.0229\\ \hline
  0.001 & 0.0007 \\
  \hline
  \end{tabular}
  \label{table7.1}
\end{table}

It is observed that stability of the method does not guarantee the accuracy. For example at $\dt = 0.01$ sec, the Leap Frog method is stable but not very accurate.

\subsection{Fully Implicit Schemes}
\label{fully_implicit_section}
\Cref{test_separate_p} is solved using three implicit schemes : The Backward Euler scheme, the Trapezoidal scheme and the Backward Differentiation Formulae (BDF2). It is observed that all the methods are unconditionally stable. For a linear system of equations (as in our case), the Trapezoidal method is equivalent to the implicit mid-point rule which is a Symplectic integrator.
At larger time steps, the results are not very accurate and the solution is damped. The following subsections give the plots of the three schemes for varying time steps.


\textbf{Fully Implicit Backward Euler Method}

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15svsBE}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE15s vs Implicit Backward Euler}
  \label{ode15svsBE}
\end{figure}

Even though the implicit Backward Euler method is unconditionally stable, it suffers from very large damping at larger time steps as observed in the \Cref{ode15svsBE}. As the time step is increased to 0.1 seconds, the solution dies out very quickly, which is unwanted for the equations represented by the Symplectic system where total energy should be conserved. The error norm is given in \Cref{table7_2}.

\textbf{Fully Implicit Trapezoidal Method}

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15s_vs_Implicit_Trap}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE15s vs Implicit Trapezoidal Method}
  \label{ode15s_vs_Implicit_Trap}
\end{figure}

As expected, the implicit Trapezoidal scheme provides better result than Backward Euler. For larger time steps (0.1 sec), the solution is not very accurate. Some phase error is observed around $x_s = 140 $, which results in large error norm for $t=0.1$ sec. The solution is more accurate for time steps around 0.01 seconds. The error norm is tabulated in \Cref{table7_2}. Impact of the phase error on the amplitude of the solution for the implicit Trapezoidal method is discussed in \Cref{Implicit_Trap_Section}


\textbf{Fully Implicit BDF2 Method}

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{BDF2vsODE15s_logcheck}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE15s vs Implicit BDF2 Method}
  \label{ode15svsBDF2}
\end{figure}

BDF2 is a multi-step second order method, and is not Symplectic in nature. As can be seen from \Cref{ode15svsBDF2}, the solution suffers from very large numerical damping at larger time steps.

In \Cref{table7_2}, the error norm for various implicit schemes at different time steps is provided. In terms of accuracy, the Trapezoidal method is better than the explicit Leap Frog scheme and BDF2 for a given time step.


\begin{table}[H]
  \centering
  \caption{Implicit Methods vs ODE15s, Error Norm}
  \vspace{0.2cm}
  \begin{tabular}{|c|c|c|c|}
  % use packages: color,colortbl
  \rowcolor{tcA}
  \hline
  Time step (sec) & Backward Euler-Norm & Trapezoidal-Norm & BDF-Norm\\  \hline
  0.1 &  0.672 & 0.470 & 0.673 \\ \hline
  0.01 & 0.672 & 0.0136 & 0.0364\\ \hline
  0.001 & 0.387 & 0.0003 & 0.007\\
  \hline
  \end{tabular}
  \label{table7_2}
\end{table}

\subsection{Semi Implicit Schemes}
\label{Semi-Implicit}

Semi Implicit schemes of the predictor corrector nature were proposed in Chapter 4 in order to make use of the existing structure of the RRB-k solver. The following algorithm is applied while using the Semi-Implicit scheme.

\begin{algorithm}                      % enter the algorithm environment
\caption{Semi Implicit Algorithm to solve \Cref{test_separate_p} }          % give the algorithm a caption
\label{alg_SI}                           % and a label for \ref{} commands later in the document
 \textbf{Input} $\vz_0$, $\vp_0$ (Initial Vector), $\dt$, itermax
\begin{algorithmic}                    % enter the algorithmic environment
    \WHILE{$t \leq tmax$}
      \WHILE{$iter  \leq itermax$}
        \IF{$iter==1$}
	  \STATE Step 1.
	  \IF {$t ==1$}
	    \STATE Advance $\vp$ explicitly using Forward Euler method on \Cref{test_separate_q} 
           \ELSE
	    \STATE Advance $\vp$ explicitly using Leap Frog method on \Cref{test_separate_q}
	   \ENDIF
        \ELSE
            \STATE Advance $\vz$ implicitly using Implicit Trapezoidal Scheme with the implicit value of $\vp$ taken from Step1.
            \STATE Correct $\vp$ implicitly using Implicit Trapezoidal Scheme with implicit value of $\vz$ taken from above Step.
        \ENDIF
        \ENDWHILE
    \ENDWHILE
\end{algorithmic}
\end{algorithm}

Here itermax represents the maximum sub-iterations performed for the corrector step of the algorithm. Apart from the time step, itermax is important for the stability considerations. Larger value of itermax should typically allow usage of large timestep. After performing numerical experiments, itermax vs maximum allowable timestep has been computed and tabulated below.

\begin{table}[H]
  \centering
  \caption{Semi Implicit Methods: itermax vs. $\dt$}
  \vspace{0.2cm}
  \begin{tabular}{|c|c|c|}
  % use packages: color,colortbl
  \rowcolor{tcA}
  \hline
  itermax & $\dt$ maximum & $\dfrac{\dt}{\dt_{cfl}}$\\ \hline
  3 & 0.019 & 1.68 \\ \hline
  5 & 0.020 & 1.77 \\ \hline
  7 & 0.021 & 1.86 \\ \hline
  9 & 0.021 & 1.86 \\ \hline
  11 & 0.021 & 1.86 \\ \hline
  19 & 0.022 & 1.95 \\ \hline
  21 & 0.022 & 1.95 \\ \hline
  23 & 0.022 & 1.95 \\ \hline
  100 & 0.022 & 1.95 \\
  \hline
  \end{tabular}
  \label{table7_3}
\end{table}

It is observed that the time step can be increased upto two times by increasing the number of iterations. 

In \Cref{ode15s_vs_SI_iter_3}, and \Cref{ode15s_vs_SI_iter_19}, the results for the semi-implicit method are plotted and compared with ODE15s.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15s_vs_SI_iter_3}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE15s vs Semi Implicit Method with itermax =3}
  \label{ode15s_vs_SI_iter_3}
\end{figure}

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{ode15s_vs_SI_iter_19}~\\[0.2cm]
  \caption{Comparison of $\zeta$ variation with length for ODE15s vs Semi Implicit Method with itermax =19}
  \label{ode15s_vs_SI_iter_19}
\end{figure}
The error norm for the Semi-implicit schemes follow the similar trend as that of the Explicit Leap Frog scheme.

% In \Cref{table7_4}, error norm for semi-implicit schemes at different time steps and itermax is provided.
% 
% \begin{table}[H]
%   \centering
%   \caption{Semi Implicit Methods vs ODE15s, Error Norm}
%   \vspace{0.2cm}
%   \begin{tabular}{|c|c|c|c|}
%   % use packages: color,colortbl
%   \rowcolor{tcA}
%   \hline
%   Time step (sec) & itermax =3 & itermax =19\\  \hline
%   0.01 &  0.0136 & 0.171 \\ \hline
%   0.019 & 0.241 & - \\ \hline
%   0.022 & - & 0.248 \\
%   \hline
%   \end{tabular}
%   \label{table7_4}
% \end{table}

% It can be seen again that increasing the number of iterations reduces the error (slightly) for a given time step. Also as compared to the explicit scheme, at the same time step, the error norm is almost equal. 
For the test problem, the parameters and the numerical techniques employed in this study, use of the Semi-Implicit methods provides little advantage over Explicit methods in terms of Stability.

\subsubsection{Impact of depth h}

The depth $h$ appears in the maximum eigenvalue of the system, and influences the stability criteria of any scheme. In order to ascertain the impact of the depth on the stability of the Semi-implicit schemes, sensitivity analysis has been carried out. The CFL number is given by \Cref{CFL_Number} and varies with the depth $h$. 

For a given itermax, it is observed that the ratio between the maximum $\dt$ achieved and the $\dt$ given by CFL number remains the same.


\section {Study of Implicit Trapezoidal Scheme}
\label{Implicit_Trap_Section}
As discussed in the previous section, the implicit Trapezoidal scheme provides a stable solution without much damping as compared to other methods. Although it suffers from the phase shift or phase error at larger time steps resulting in larger error norms.


In this section, the impact of phase difference on the amplitude of the solution for the coupled set of equations and non-coupled set of equations is discussed.


\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_NoCoupling_t=20}~\\[0.2cm]
  \caption{Implicit Trapezoidal method without Coupling at t = 20sec for various time steps}
  \label{Implicit_Trap_NoCoupling_t=20}
\end{figure}


In \Cref{Implicit_Trap_NoCoupling_t=20} and \Cref{Implicit_Trap_NoCoupling_t=50}, the values of variable $\zeta$ are plotted at two differnt time intervals of 20sec and 50 sec respectively. The system of equations have been decoupled by setting $h=0$ and $g=0$. No phase error is observed for this case even when large time steps are used. As the initial solution do not contain sharp discontinuities, Central discretization method does not introduce  oscillations and the phase errors are minimal. This changes when we introduce the coupling between the two equations by setting $h=50$ and $g=9.81$.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_NoCoupling_t=50}~\\[0.2cm]
  \caption{Implicit Trapezoidal method without Coupling at t = 50sec for various time steps}
  \label{Implicit_Trap_NoCoupling_t=50}
\end{figure}

For the coupled equations (\Cref{Implicit_Trap_Coupling_t=22_24_26}), it is observed that the amplitude of the solution does not remain constant due to the impact of the source term, and varies with the phase. This impact can be seen in \Cref{Implicit_Trap_Coupling_t=22_24_26} where the benchmark solution is plotted at three different time intervals. Without any coupling between the equations, the solution should travel with the wave speed without changing the shape or the amplitude. With the coupling, a change in amplitude is seen at different time intervals.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_Coupling_t=22,24,26}~\\[0.2cm]
  \caption{Solution of coupled equation at t = 20,22 and 24 sec with Matlab ODE integrator}
  \label{Implicit_Trap_Coupling_t=22_24_26}
\end{figure}

This change in the amplitude with the phase for the coupled system causes large error norms for Implicit Trapezoidal method. As Implicit Trapezoidal method induces a phase error, the change in the phase results in a change in amplitude. This is reflected in the \Cref{Implicit_Trap_Coupling_t=20} 

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_Coupling_t=20}~\\[0.2cm]
  \caption{Implicit Trapezoidal method with Coupling at t = 20 sec for various time steps}
  \label{Implicit_Trap_Coupling_t=20}
\end{figure}


For time step $\dt =0.01 s$, the solution is accurate and does not suffer from phase error. For larger timer steps, both phase change and amplitude variation is observed. It is important to observe that, for $\dt=0.1 s$, the solution appears to be already damped at time 20s. If there was damping in the solution, then it would be greater for larger time steps, which is not seen for $\dt=1 s$. Instead the amplitude at $\dt=1 s$ is greater than the benchmark solution.

It could be argued that the larger $\dt$ is impacting the stability of the solution. In order to check this, the solution is plotted for time 100 sec in \Cref{Implicit_Trap_Coupling_t=100}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_Coupling_t=100}~\\[0.2cm]
  \caption{Implicit Trapezoidal method with Coupling at t = 100 sec for various time steps}
  \label{Implicit_Trap_Coupling_t=100}
\end{figure}

It is observed that the solution remains stable, and only phase error and amplitude difference is observed for larger time steps.

\subsection{Higher Order Schemes and Impact on Phase Error}

Higher order (Two stage, Fourth ordered) Symplectic Implicit Runge Kutta scheme has been implemented to analyze the impact on the phase error due to the larger time steps. These methods are computationally more intense, but provides better accuracy.
Given the coupled set of equations as:

\begin{equation*}
 \dfrac{d \textbf{q}}{dt} + L\textbf{q} =0 
\end{equation*}

we have

\begin{equation}
\label{RK4}
 \dfrac{q^{n+1}-q^n}{\dt} = \dfrac{1}{2} (K_1 + K_2)
\end{equation}

where 
\begin{align*}
 K_1 &= -L (q^n + \dt/4 K_1 + (1/4 + \sqrt{3}/6) \dt K_2) \\
 K_2 &= -L (q^n + (1/4 - \sqrt{3}/6)\dt K_1 + 1/4\dt K_2) \\
\end{align*}

At every time step, above system of equation is simultaneously solved and then plugged into \Cref{RK4}. Please note that the coefficients are chosen such that the method remains Symplectic in nature.

The results are shown in \Cref{RK4_image}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{RK4}~\\[0.2cm]
  \caption{Higher order Runge Kutta method at time 100s}
  \label{RK4_image}
\end{figure}

It is observed that time step $\dt =0.1s$, the solution is accurate and there is no phase error. As compared to the Implicit Trapezoidal method, higher order Runge Kutta method is more accurate. Phase error is observed when the time step is further increased. 

\subsection{Impact of Initial Conditions}

In a real case scenario, the incoming waves are generated by mixing various waves of different frequencies and amplitudes. The initial condition in this case is chosen as:

\begin{equation}
 \zeta^0  (x) = 0.25( \sum_{i=0}^{3} cos(\dfrac{2^i4\pi x}{L} )  
\end{equation}

For the initial conditions with single wave, the solution was accurate for $\dt =0.01s$. The solution for same $\dt =0.01$ but the mixed initial conditions is plotted in \Cref{Implicit_Trap_Coupling_MixWaves_t=100_dt=0.01}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_Coupling_MixWaves_t=100_dt=001}~\\[0.2cm]
  \caption{Implicit Trapezoidal Method with mixed initial conditions, $\dt$ =0.01s}
  \label{Implicit_Trap_Coupling_MixWaves_t=100_dt=0.01}
\end{figure}

It is observed that the solution is not accurate anymore. This can be attributed to high frequency waves. The phase error for Implicit Trapezoidal schemes is generally proportional to square of the frequency of the wave.

Instead of a mix of waves, simulation is carried with initial condition as $\zeta^0  (x) = cos(\dfrac{32\pi x}{L} ) $. The result for $\dt =0.01s$ is plotted in \Cref{Implicit_Trap_HighFrequency}. Phase error is observed even for a single initial wave with high frequency at $\dt=0.01s$. This implies that for high frequency waves, a smaller time step is required to maintain accuracy. 

Spatial discretization also influences the solution with the high frequency waves. As a rule of thumb, a minimum of 15 grid points per wavelength should be used to accurately represent the wave. In our study, we have used 12.5 gridpoints per wavelength. Increasing the number of grid points reduces $\Delta x$, and thus will result in lower CFL limit. In general, it is the low frequency waves, with larges wavelengths which influences the Ships most. The $\Delta x$ should be chosen to represent these low frequency waves accurately. 

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_HighFrequency}~\\[0.2cm]
  \caption{Implicit Trapezoidal Method with high frequency inlet conditions, $\dt$ =0.01s}
  \label{Implicit_Trap_HighFrequency}
\end{figure}


\subsection{Impact of Boundary Conditions}


Periodic boundary conditions are prone to accumulation of phase error as the time grows. In this section, we try to analyze the impact on the solution when different boundary conditions are used. 

We use the following:

\begin{itemize}
 \item $\zeta(0,t) = 1.0$ : Dirichlet Boundary Condition.
 \item For numerics, $\zeta_{N+1} = \zeta_{N}$ where N is the node number at right boundary. This is also equivalent to the Neumann boundary condition.
 \item Neumann Boundary Condition for $\varphi$
\end{itemize}

The results are shown in \Cref{Implicit_Trap_BC}. Again, phase errors are observed in the case of larger time steps.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_BC}~\\[0.2cm]
  \caption{Impact of boundary condition on Implicit Trapezoidal Method}
  \label{Implicit_Trap_BC}
\end{figure}


\newpage
\subsection{Summary}

A summary of the methods is provided below.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Comparison}~\\[0.2cm]
  \caption{Comparison of various time integration schemes and stability estimates}
  \label{comparison}
\end{figure}
\newpage
\section{Impact of Spatial Discretization}

As discussed in \Cref{benchmark}, wiggles/oscillations are observed in $\zeta$ when the initial condition has a discontinuity. Following are the possible reasons.
\begin{itemize}
 \item Spatial Discretization
  \item Coupling of equations
\end{itemize}

In order to analyze the impact of the coupling on the oscillations, \Cref{test_separate_p} is decoupled by setting $h=0$ and $g=0$. This results into two decoupled wave equations, whose analytical solution is available. Let the initial condition be given as $\zeta^0(x)$ and $\varphi^0(x)$. Then the solution at $(x,t)$ is given by :

\begin{align*}
 \zeta (x,t) & = \zeta^0(x-Ut) \\
 \varphi (x,t) & = \varphi^0(x-Ut)
\end{align*}

Based on the above equations, if the initial condition has a discontinuity, then it should propagate with time without changing any shape. Spatial discretization usually impacts the solution of the wave equation in the case of discontinuities or shock-waves. 

Again, equations are discretized using Central discretization and solved using the Matlab benchmark code (ODE15s). The results obtained for $\zeta$ are given in \Cref{Central_without_h_g}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Central_without_h_g}~\\[0.2cm]
  \caption{$\zeta$ for decoupled wave equation after one time period with Central discretization and with discontinuity in the initial condition}
  \label{Central_without_h_g}
\end{figure}

Wiggles and oscillations are still present in the solution, which leads us to explore the other possibility of Spatial discretization. Whereas Central discretization causes oscillations around the discontinuity, error introduced by the Upwinding discretization be classified into following:

\begin{itemize}
 \item Smearing: The length of Smearing after the discontinuity.
 \item Phase Error
\end{itemize}

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Implicit_Trap_Upwinding_Error}~\\[0.2cm]
  \caption{Classification of Errors induced by Upwinding Discretization}
  \label{Implicit_Trap_Upwinding_Error}
\end{figure}

In \Cref{Implicit_Trap_Upwinding_Error}, the anlaytical solution is shown as the black line, whereas the solution from the Upwinding discretization is shown in red.

Instead of Central discretization, Upwinding scheme is now used to discretize \Cref{Reference_Test_case}. This results in following discretization stencils:

\begin{align*}
 \Szz  &: \begin{bmatrix} \dfrac{-U}{\Delta x} &&  \dfrac{U}{\Delta x} && 0 \end{bmatrix}  \\  
 \Szp  &: \begin{bmatrix}  \dfrac{h} {\Delta x^2} && \dfrac{-2 h} {\Delta x^2} && \dfrac{h} {\Delta x^2} \end{bmatrix} \\
 \Spz  &: \begin{bmatrix} 0 && g && 0 \end{bmatrix} \\
 \Spp  &: \begin{bmatrix} \dfrac{-U}{\Delta x} &&  \dfrac{U}{\Delta x} && 0 \end{bmatrix} \\  
\end{align*}

Using Matlab ODE15s and the above discretization stencils, the solution is obtained is plotted in \Cref{Upwinding_without_h_g}:

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Upwinding_without_h_g}~\\[0.2cm]
  \caption{$\zeta$ for decoupled wave equation after one time period with Upwinding discretization}
  \label{Upwinding_without_h_g}
\end{figure}

In \Cref{Upwinding_without_h_g}, we observe that the solution is damped, but no oscillations or wiggles are present in the solution. Also, larger damping at $x=0$ is observed which can be attributed to the discontinuity.

Damping is a characteristic of the Upwinding discretization and can be removed by adopting higher order or higher resolution methods. In order to confirm that the wiggles are generated only because of the central discretization scheme and not because of the coupling, we run the simulation with the Upwinding discretization and by setting $h=50$ and $g=9.81$, which again gives us coupled set of equations.
Again running the simulation with Matlab ODE15s integrator, we obtain the following result:

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{Upwinding_with_h_g}~\\[0.2cm]
  \caption{$\zeta$ for Coupled wave equation after one time period with Upwinding discretization}
  \label{Upwinding_with_h_g}
\end{figure}

In \Cref{Upwinding_with_h_g}, for the coupled equations, the Upwinding discretization removes the wiggles and also results in a smooth but damped solution.

We conclude that the wiggles are generated in the simulation if a discontinuity is present in the initial condition, with the Central discretizations scheme. However, it does not impact the stability of the Semi-implicit schemes.

\newpage
\subsection{Impact of Grid size on Accuracy}

One of the goals of the current project is to make the mesh finer near the regions surrounding a ship or obstacles. As we have seen that varying the time step impacts the accuracy of the solution, we will try to assess the time step requirements for similar accuracy for different grid sizes.

Simulations are run with the implicit Trapezoidal scheme, and varying grid spacing and adjusting time step. As the Implicit Trapezoidal method is second order accurate, reduction in grid size implies reduction in time step if the CFL number is not changed. Thus the accuracy of the solution improves by $\dt^2$, and a better benchmark is required to compare the solution. Again the Matlab ODE15s method is used for the benchmarking along with an absolute error estimate of $1e^{-9}$ as compared to the previous absolute error of $1e^{-6}$. The results obtained are indicated in \Cref{X_comparison}

\begin{table}[htbp]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{$\dx$} & \textbf{CFL Number} & \textbf{Error Norm} \\ \hline
\multicolumn{ 1}{|c|}{2} & 1.8 & 0.717622 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.9 & 0.235133 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.4 & 0.053416 \\ \hline
\multicolumn{ 1}{|c|}{1} & 1.8 & 0.201493 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.9 & 0.046899 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.4 & 0.015783 \\ \hline
\multicolumn{ 1}{|c|}{0.5} & 1.8 & 0.011764 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.9 & 0.0063474 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.4 & 0.004944 \\ \hline
\multicolumn{ 1}{|c|}{0.25} & 1.8 & 0.005404 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.9 & 0.005445 \\ \cline{ 2- 3}
\multicolumn{ 1}{|c|}{} & 0.4 & 0.003001 \\ \hline
\end{tabular}
\end{center}
\caption{Error Norm for various Grid sizes and CFL numbers}
\label{X_comparison}
\end{table}


For the same CFL number, an error reduction from 0.71 to 0.2 is observed when the grid spacing is reduced from 2 to 1. This approximate factor of 4 is attributed to the reduction in time step by two times and the second order temporal accuracy. Another observation is that, for the combination of small grid size ($\Delta x = 0.25$) and small $\dt$, the error norm does not improve with further reduction in $\dt$.

\section{Stability Analysis with Vertical Structure $\psi$ included}

Based on the observations in \Cref{fully_implicit_section}, the fully Implicit Trapezoidal scheme provides the best stability properties for the test case given by \Cref{Reference_Test_case}. We still need to evaluate the impact of the vertical structure $\psi$ which is given by an elliptic equation and does not include a time derivative term. To assess the impact of this additional variable and equation, the following test case has been considered.

\begin{subequations}
\label{Reference_Test_case_2}
\begin{equation}
\frac{\partial \zeta}{\partial t} +\textbf{U} \dfrac{\partial \zeta}{\partial x} + h \dfrac{\partial^2 \varphi}{\partial x^2} -  h\mathcal{D} \dfrac{\partial^2 \psi}{\partial x^2}= 0,
\end{equation}
\begin{equation}
\frac{\partial \varphi}{\partial t} + \textbf{U}\dfrac{\partial \varphi}{\partial x} + g\zeta = 0,
\end{equation}
\begin{equation}
\mathcal{M}\psi +  h D \dfrac{\partial^2 \varphi}{\partial x^2} -  \NN \dfrac{\partial^2 \psi}{\partial x^2} = 0,
\end{equation}
\end{subequations}

Three additional discretization matrices are formulated. The model parameters $\mathcal{M} , \NN \mathcal{D}$ are derived from the vertical structure shape. For the linearized parabolic model we obtain, 
\begin{equation}
 \mathcal{D} = \dfrac{1}{3}h , \quad  \mathcal{M} = \dfrac{1}{3}h, \quad \textbf{and } \NN = \dfrac{2}{15}h^3
\end{equation}

The discretization matrices are given as:

\begin{align*}
 \Szz  &: \begin{bmatrix} \dfrac{-U}{2 \dx} && 0  && \dfrac{U}{2 \dx}\end{bmatrix}\\
 \Szp  &: \begin{bmatrix}  \dfrac{h} {\dx^2} && \dfrac{-2 h} {\dx^2} && \dfrac{h} {\dx^2} \end{bmatrix} \\
 \Spz  &: \begin{bmatrix} 0 && g && 0 \end{bmatrix} \\
 \Spp  &: \begin{bmatrix} \dfrac{-U}{2 \dx} && 0  && \dfrac{U}{2 \dx}\end{bmatrix} \\
 \Szs  &: \begin{bmatrix}  \dfrac{h^2} {3\dx^2} && \dfrac{-2 h^2} {3\dx^2} && \dfrac{h^2} {3\dx^2} \end{bmatrix} \\
 \Ssp  &: \begin{bmatrix}  \dfrac{h^2} {3\dx^2} && \dfrac{-2 h^2} {3\dx^2} && \dfrac{h^2} {3\dx^2} \end{bmatrix} \\
 \Sss  &: \begin{bmatrix}  \dfrac{-2h^3} {15\dx^2} && \dfrac{4h^3} {15\dx^2} && \dfrac{-2h^3} {15\dx^2}\end{bmatrix} \\
\end{align*}

Periodic boundary conditions for $\psi$ has been assumed. The model parameters $h$, $U$ and $g$ remain the same as in the case of two equation model.

\subsection{Benchmark}
\label{benchmark1}
Equations for $\zeta$ and $\varphi$ are coupled and represented as a Lumped model. It contains terms from $\psi$ which is explicitly not known as it depends on $\varphi$.

Let the lumped equation be represented as :

\begin{equation}
 \dfrac{d\vec q}{dt} + L \vec q = C \vec \psi
\end{equation}

where $L=\begin{bmatrix} \Szz && \Szp \\ \Spz && \Spp \end{bmatrix}$ and $C=-\begin{bmatrix} \Szs \\ \Sps \end{bmatrix}$


The set of differential algebraic equations given by \Cref{Reference_Test_case_2} can be solved by the Matlab ODE15s solver, an implicit method capable of solving stiff problems and DAE's. The derivatives in time are described as :

\begin{equation}
 \begin{bmatrix} \dfrac{d\vec q}{dt} \\ \dfrac{d \vs}{dt} \end{bmatrix} = \begin{bmatrix} -L \vec q + C \vs \\ - \Ssp \vp - \Sss \vs \end{bmatrix} 
\end{equation}

In order to solve the DAE by Matlab ODE15s, an additional input given by $M = \begin{bmatrix} I_{2x_s \times 2x_s} && 0 \\ 0 && 0\end{bmatrix}_{3x_s \times 3x_s}$ is required, where $x_s$ is the number of grid points. This matrix M, also called as the Singular Matrix characterizes the 
Also, one has to be careful while chosing the initial conditions for the system of DAE's. An unique solution is obtained when the initial conditions also satisfy the the algebraic equation. We have chosen the initial values of $\vs$ and $\vp$, such that they satisfy the algebraid equation.

% Following solution is obtained after running ODE15s for 100s.
% 
% \begin{figure}[H]
% \centering
% \includegraphics[width=0.8\textwidth]{ode15s_DAE}~\\[0.25cm]
% \caption{$\zeta$ variation with length at different times using ODE15s for the DAE}
% \label{ode15s_DAE}
% \end{figure}

\subsection{Chorin's Projection technique}

In order to solve this coupled equation, a projection method similar to Chorin's projection method \cite{Chorin} has been used. Chorin's Projection scheme was originally used to solve the incompressible unsteady state Navier Stokes equations. The key advantage of the projection method is that the computation of the dependent variables are decoupled (Pressure and velocity in the case of Navier Stokes equations). Typically the algorithm consists of two stage fractional step schene, a method which uses multiple calculation steps for each numerical time-step. 

The steps are split as following:

\begin{itemize}
 \item Ignore the $\psi$ term from the lumped equation, and solve the system $\dfrac{d\vec q}{dt} = -L \vec q $ with the Implicit Trapezoidal method. Denote the solution after one time step as $\vec q*$
 \item Compute $\vs_{n+1}$ from $\Ssp \vp + \Sss \vs =0$ by solving the system of equations where $\vp$ is derived from $q*$.
 \item Compute $\vec q_{n+1}$ as $\vec q_{n+1} = \vec q* + \dt C \vs_{n+1}$
\end{itemize}

In \Cref{ODE15s_vs_Chorin_Check}, the solution obtained from Chorin's projection scheme for different time steps is compared with the Benchmark solution.

\begin{figure}[H]
  \centering
  \includegraphics[width=0.8\textwidth]{ODE15s_vs_Chorin_Check}~\\[0.25cm]
  \caption{$\zeta$ variation with length at different times using ODE15s using Chorin's Projection}
  \label{ODE15s_vs_Chorin_Check}
\end{figure}

% In order to compare the solutions at different time intervals, the error norm with respect to the benchmark ODE15s solution is tabulated.
% 
% \begin{table}[H]
%   \centering
%   \caption{Chorin Projection vs ODE15s, Error Norm}
%   \vspace{0.2cm}
%   \begin{tabular}{|c|c|}
%   \hline
%   \rowcolor{tcA}
%   Time step (sec) & Error Norm\\
%   \hline
%   0.01 & 0.175 \\ \hline
%   0.1 & 0.543 \\
%   \hline
%   \end{tabular}
%   \label{table7_6}
% \end{table}

As the time step is reduced, the solution is more close to the Benchmark solution. As the time step is increased, a larger phase error can be observed. High error norms for the solution at the larger time steps can be attributed to this phase error. 

% In order to check the Symplectic behaviour of the solution, we compute the total strenth of the Hamiltonian which is described in Chapter 2 in detail. This is done by simply adding the values of $\zeta$ at all grid points. Following is the comparison:
% 
% \begin{table}[H]
%   \centering
%   \caption{Chorin Projection vs ODE15s, Hamiltonian Strength}
%   \vspace{0.2cm}
%   \begin{tabular}{|c|c|}
%   \hline
%   \rowcolor{tcA}
%   Case & Hamiltonian Stength\\
%   \hline
%   Initial condition  & -3.5577 \\ \hline
%   ODE15s benchmark & -3.5577 \\ \hline
%   Chorin Projection $\dt = 0.01s$ &  -3.5577 \\ \hline
%   Chorin Projection $\dt = 0.1s$ &   -3.5577 \\
%   \hline
%   \end{tabular}
%   \label{table7_7}
% \end{table}
% 
% Results in \Cref{table7_7} prove that the implicit trapezoidal method conserves the energy of the system, even at higher time steps where large oscillations and phase shift are observed. This is an important conclusion, as the system under consideration is a Hamiltonian system, and Symplectic nature of the numerical scheme is a requirement. Also, with the Chorin's projection scheme, we obtain an unconditionally stable scheme. Similar results are observed for time step $\dt =1 s$.

\section{Stability Analysis for Two Dimensional Case}

In order to analyze the stability of the complete system, it is required to study the two dimensional case which results in a Pentadiagonal matrix instead of a Tridiagonal matrix.

Following system of equations are considered as the test case:

\begin{subequations}
\label{Reference_Test_case_3}
\begin{equation}
\frac{\partial \zeta}{\partial t} +\textbf{U}_x \dfrac{\partial \zeta}{\partial x} + \textbf{U}_y \dfrac{\partial \zeta}{\partial y}  + h (\dfrac{\partial^2 \varphi}{\partial x^2} + \dfrac{\partial^2 \varphi}{\partial y^2})-  h \mathcal{D} ( \dfrac{\partial^2 \psi}{\partial x^2} - \dfrac{\partial^2 \psi}{\partial y^2})= 0,
\end{equation}
\begin{equation}
\frac{\partial \varphi}{\partial t} + \textbf{U}_x\dfrac{\partial \varphi}{\partial x} + \textbf{U}_y\dfrac{\partial \varphi}{\partial y} + g\zeta = 0,
\end{equation}
\begin{equation}
\mathcal{M}\psi +  h D ( \dfrac{\partial^2 \varphi}{\partial x^2} +\dfrac{\partial^2 \varphi}{\partial y^2}) -  \NN (\dfrac{\partial^2 \psi}{\partial x^2} + \dfrac{\partial^2 \psi}{\partial y^2}) = 0,
\end{equation}
\end{subequations}

Finite Volume method based on Central discretization has been used for the spatial discretization of the system of equations. A rectangular grid has been assumed with possibility of different grid size for the $X$ and $Y$ direction. The discretization matrices are Pentadiagonal matrices represented by 5-point stencil as described in Chapter 2 in detail.

It is important to consider the boundary conditions while constructing the Pentadiagonal matrices. In our case, periodic boundary conditions in both X and Y directions have been used. 

The model parameters are given below for the Two dimensional case:

\begin{itemize}
  \item $L_x$ =100 m
  \item $L_y$ =100 m
  \item $N_x$ = 20 	Number of grid points
  \item $N_y$ = 20 	Number of grid points
  \item $U_x$ = $U_y$ = 1.0 m/s
  \item h = 50 m
  \item time-max = 50 sec
  \item Initial Condition: $ \zeta^0  (x,y)$ = cos($\dfrac{4\pi x}{L_x}$)  
  \item $ \varphi^0  (0,x,y) = 1$
  \item $ \psi^0(x,y) =0$
\end{itemize}

The initial condition for $\zeta$ is shown below:

\begin{figure}[H]
  \centering
  \includegraphics[width=0.8\textwidth]{initial_2D}~\\[0.25cm]
  \caption{Initial conditions for $\zeta$ for the two Dimensional case}
  \label{initial_2D}
\end{figure}


It can be seen in \Cref{initial_2D} that the initial conditions are periodic in nature but as the grid is very coarse, the initial condition is not very smooth in nature. The way currently the periodic boundary conditions are implemented, a non-smooth initial condition impacts the periodicity of the solution with time if we use the Central Discretization scheme.

The equations in lumped form for the benchmark solution and the Chorin's projection remains the same. Instead of the Tridiagonal matrices, we have pentadiagonal matrices and block matrices. The grid is ordered lexicographically in Column major format (Similar to Fortran ordering).

Stability, accuracy and the Symplectic nature of the solution from Chorin's projection method is compared with the Benchmark solution. The Implicit Trapezoidal method has been used as the substep of the Chorin's projection method. Please note that, as the mesh size has been kept large to keep the problem size small, the level of accuracy expected is low for moderate time steps (~0.1 sec).

Following  figures show the comparison of the Chorin's projection method vs. the benchmark solution for two different time steps.

\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{twod_dt001}~\\[0.25cm]
  \caption{$\zeta$ variation using ODE15s vs. Chorin's Projection for $\dt=0.01s$}
  \label{twod_dt0.01}
\end{figure}


\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{twod_dt01}~\\[0.25cm]
  \caption{$\zeta$ variation using ODE15s vs. Chorin's Projection for $\dt=0.1s$}
  \label{twod_dt0.1}
\end{figure}

It is observed that the solution is accurate for the smaller time steps, and phase error can be seen for $\dt=0.1 s$. Also, the benchmark solution is not periodic in nature. The implementation of the periodic boundary condition results in such discrepancy. In the current implementation, the value at Node 0 is set equal to Node $N$. If it is changed from Node $N$ to Node $N-1$, then the solution becomes periodic. Also, if simulations are carried with finer grid, the periodic behaviour is seen again as the value at Node $N$ and Node $N-1$ are no longer different. The solution is shown in \Cref{twod_dt1}.


\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{twod_dt001_n_50}~\\[0.25cm]
  \caption{$\zeta$ variation using ODE15s vs. Chorin's Projection for $\dt=0.01s$ and finer grid}
  \label{twod_dt001_n_50}
\end{figure}



\begin{figure}[H]
  \centering
  \includegraphics[width=1.0\textwidth]{twod_dt1}~\\[0.25cm]
  \caption{$\zeta$ variation using ODE15s vs. Chorin's Projection for $\dt=1s$ and finer grid}
  \label{twod_dt1}
\end{figure}

The system is stable for larger time steps, which gives us confidence in proceeding with Chorin's Projection scheme for the real case.
% 
% \Cref{table7_8} comtains the errors-norm associated for different time steps. Clearly at larger time steps, due to the phase shift, error norm is not the correct way of analyzing results.
% 
% \begin{table}[H]
%   \centering
%   \caption{Chorin Projection vs ODE15s- Two dimensional, Error Norm}
%   \vspace{0.2cm}
%   \begin{tabular}{|c|c|}
%   \hline
%   \rowcolor{tcA}
%   Time step (sec) & Error Norm\\
%   \hline
%   0.01 & 0.0254 \\ \hline
%   0.1 & 0.931 \\ \hline
%   1   &	 1.04 \\
%   \hline
%   \end{tabular}
%   \label{table7_8}
% \end{table}
% 
% \textbf{Symplectic Nature}
% 
% In order to check the symplectic nature of the solution, the values of $\zeta $ at all grid points are simply added for the various solutions and the initial conditions. As we have used the implicit Trapezoidal method with the Chorin's projection, we obtain the constant energy of -34.4636 for all the solutions.
% 
